Numerical study of large-scale vorticity generation in shear-flow turbulence 
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Simulations of stochastically forced shear-flow turbulence in a shearing-periodic domain are used to study 
the spontaneous generation of large-scale flow patterns in the direction perpendicular to the plane of the shear. 
Based on an analysis of the resulting large-scale velocity correlations it is argued that the mechanism behind 
this phenomenon could be the mean-vorticity dynamo effect pioneered by Elperin, Kleeorin, and Rogachevskii 
in 2003 (Phys. Rev. E 68, 01631 1). This effect is based on the anisotropy of the eddy viscosity tensor. One of 
its components may be able to replenish cross-stream mean flows by acting upon the streamwise component of 
the mean flow. Shear, in turn, closes the loop by acting upon the cross-stream mean flow to produce stronger 
streamwise mean flows. The diagonal component of the eddy viscosity is found to be of the order of the rms 
turbulent velocity divided by the wavenumber of the energy-carrying eddies. 

PACS numbers: PACS Numbers : 47.27.tb, 47.27.ek, 95.30.Lz 



I. INTRODUCTION 



The imperfect analogy between the induction equation and 
the vorticity equation has always raised questions regarding 
the extent of this analogy. While it is well-known that the av- 
eraged induction equation for the mean magnetic field admits 
self-excited solutions for a turbulent flow with helicity, anal- 
ogous solutions to the averaged vorticity equation only exist 
in the compressible case |T, z]. An exception is the case of 
flows that are driven by a non-Galilean invariant forcing func- 
tion, which can give rise to the so-called anisotropic kinetic 
alpha effect IHUIllI^. This effect produces mean flows that 
are helical and of Beltrami type. Another example of mean 
flow generation is the A effect [7, 8], whereby large-scale non- 
uniform flows can be produced in rotating anisotropic turbu- 
lence. 

In the last few years another example has emerged, where 
the analogy between vorticity and induction equations is more 
striking. This example applies to the case of shear-flow tur- 
bulence. In fact, it has been argued that large-scale magnetic 
field generation is possible via the shear-current effect that re- 
sults from non-vanishing off-diagonal components of the tur- 
bulent magnetic diffusivity tensor |9, 10]. This effect predicts 
large-scale field generation in homogeneous shear-flow turbu- 
lence with non-helical driving, which has indeed been seen in 
several simulations 11311 . However, there is the prob- 

lem that, according to the test-field method, the sign of the 
relevant component of the turbulent magnetic diffusivity ten- 
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sor was found to be incompatible with that required for the 
shear-current dynamo 1 12]. On the other hand, the analogous 
hydrodynamic effect has not yet been explored in sufficient 
detail. This effect may explain the generation of large-scale 
vorticity in homogeneous shear-flow turbulence and was first 
studied analyticall y iii a seminal paper by Elperin, Kleeorin, 
and Rogachevskii 01411 . Several recent studies discuss numer- 
ical evidence for the spontaneous formation of mean vorticity 
[11, 12, 13]. In those papers the main objective is to study 
the generation of large-scale magnetic fields by shear-flow tur- 
bulence, while the simultaneous generation of mean vorticity 
was merely an additional (but interesting) complication. On 
the other hand, in view of the disappointing experience when 
trying to verify the operation of the shear-cuiTent dynamo us- 
ing the test-field method, one should be careful in view of 
earlier negative results concerning both the shear-current 
effect and the mean-vorticity dynamo effect. The aim of this 
paper is therefore to discuss turbulent shear flow simulations 
without magnetic fields in order to demonstrate the existence 
of the mean-vorticity dynamo and to analyze its connection 
with the eddy viscosity tensor in more detail. 

Following earlier work lfTll[T2ilT3^ . periodic boundary con- 
ditions are used in the streamwise direction and in the direc- 
tion perpendicular to the plane of the shear, while shearing- 
periodic boundary conditions are used in the cross-stream di- 
rection. This means that mass and mean momentum are con- 
served. Furthermore, if a large-scale flow emerges, it will also 
be periodic coiTesponding to a simple sine wave. The mean 
vorticity is therefore also a long-wavelength sine wave. How- 
ever, although the original analysis was based on mean vor- 
ticity, we discuss in the following mainly the mean velocity, 
because the corresponding equations are simpler and more in- 
tuitive. 

For a proper analysis of the hydrodynamic mean-vorticity 
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dynamo effect one would need to proceed analogously to the 
hydromagnetic case where it was possible to determine all rel- 
evant components of the turbulent magnetic diffusivity ten- 
sor using the test-field method. One would then need to de- 
termine all relevant components of the eddy viscosity tensor. 
However, in the absence of a properly developed "test-flow" 
method for hydrodynamics, we must resort to more primitive 
measures for estimating components of the eddy viscosity ten- 
sor. Using decay calculations of a large-scale velocity struc- 
ture, it was found that eddy viscosity, v^, and turbulent mag- 
netic diffusivity, -q^, are approximately equal, i.e. w T^t, and 
around (0.8. ..0.9) x Wrms/fcf 116]. Here, fcf is the wavenum- 
ber corresponding to the scale of the energy-carrying eddies 
and Wrms is the rms velocity of the turbulence. On the other 
hand, a more accurate determination of rjt led recently to the 
'Ht = Vto = ""rms/ (3fct), where ?7to is just a reference value. 
In this paper we use an analogously defined reference value, 
i^to = Urms/(3fcf), but note that there is no strong case for 
assuming that will be close to vtQ. 



II. THE MODEL 

In the present work we consider weakly compressible sub- 
sonic turbulence in the presence of a linear shear flow. 



U" 



(0,^x,0), 



(1) 



so X is the cross-stream direction, y is the streamwise direc- 
tion, and z is the direction perpendicular to the plane of the 
shear flow. Since the effect of temperature changes is not im- 
portant in this context, we consider an isothermal equation of 
state. In the following we work with the departures from this 

mean flow, so the total velocity is J7 +U, and the governing 
equations for U are then Il2ll 



Dlnp 



-V-t/, 



(2) 



(3) 



where T) /T)t — d / dt + {U + [/) • V is the advective deriva- 
tive with respect to the total velocity, Cs is the isothermal 
sound speed, here considered as constant, p is the mass den- 
sity, / is a random forcing function, -Fvisc = p^^V • 2pvS 
is the viscous force, and Sy = \{Ui,j + Uj^i) — ^Sij'V ■ U 
is the traceless rate of strain tensor and commas denote partial 
derivatives. 

The forcing function is (5-correlated in time and consists 
of random plane waves with wavevectors k in the interval 
4.5 < k/ki < 5.5 [17]. During each time step, / is a sin- 
gle transverse (solenoidal) plane wave proportional to fc x e, 
where the wavevector k is taken randomly from a set of pre- 
defined vectors with components that are integer multiples of 
2tt/L and whose moduli are in a certain interval around an av- 
erage value, ( I A; I ) , which we denote by fcf , and e is an arbitrary 
random unit vector not aligned with k. The corresponding 



scale, 27r/fcf, is referred to as the energy-carrying scale of the 
turbulence. Moreover, the time dependence of / is designed 
to mimic ^-correlation, which is a simple and commonly used 
form of random driving [T^] . 

There are two important dimensionless control parameters, 
the Reynolds number Re and the shear parameter Sh, 



Re = Urais/{l^k{), Sh = 5/(Urmsfct), 



(4) 



that quantify the intensity of turbulence and shear, respec- 
tively. We note that the values of Re and Sh cannot be chosen 
a priori due to the strong effect that the vorticity dynamo has 
on the value of itrms in the saturated state. Thus we always 
refer to values of Uims, Re ,and Sh that apply to the situation 
where the vorticity dynamo is absent, i.e. early stages of the 
run or a non-shearing simulation. The ratio of the size of the 
domain, L, to the size of the energy-carrying scale is also an 
important control parameter that we call the scale separation 
ratio, written here as kfL /2tt = kf/ki, where fci — 2tt/L is 
the smallest wavenumber that fits into the domain. 

We employ the PENCIL Code |18] with sixth-order finite 
differences in space and a third order time stepping scheme. 
We use triply-periodic boundary conditions, except that the x 
direction is shearing-periodic, i.e. 



U{-^Lx,y,z,t) = U{^Lx,y + LxSt,z,t). 



(5) 



This condition is routinely used in numerical studies of shear 
flows in Cartesian geometry lfl9[l20ll . 



III. RESULTS 

The initial velocity is zero, but the volume forcing drives a 
random flow that soon develops a turbulent cascade where the 
spectral energy follows an approximate k^^^^ inertial range 
between the forcing wavenumber fcf and some dissipation 
wavenumber fcd — {S^ /v^)^/^. 

In Fig.[T]we show images of the streamwise component of 
U at the periphery of the computational domain from a run 
with Re K, 100 and Sh « -0.2 (hereafter Run A). At early 
times the velocity pattern is dominated by structures whose 
scale is comparable with the forcing scale, which is about one 
fifth of the domain size. However, at later times there is a ten- 
dency to produce large-scale flow patterns with a long wave- 
length variation in the z direction. This flow pattern tends to 
be unstable and keeps disappearing and reappearing. This is 
seen also for other runs with smaller Reynolds number. 

Given the systematic variation in the z direction, it is useful 
to consider averages over the x and y directions, denoted in 
the following by overbars. So, U = U{z, t) depends only on 
z and t. Figure |2] shows Ux and Uy as functions of time and 
z. In Fig. [3] we plot the z dependence of Ux and Uy at a time 
near the maximum vorticity. Note that the amplitude of Uy 
is about 4 times as big as that of Ux, and that the two fields 
are essentially in phase. The fact that Ux and Uy are in phase 
is an immediate consequence of the fact that 5* < 0, and that 
there is a minus sign in front of 5* in Eq. (|2|i. 
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FIG. 1: (Color online) Representation of Uy on the periphery of the computational domain for Run A at six different times showing the 
occasional generation of large scale flow patterns with a systematic variation in the z direction. Dark (blue) shades refer to negative values of 
Uy while light (yellow) shades refer to positive values. Note that at time tCski = 900 (corresponding to turmskf ~ 480) the orientation of the 
flow pattern in the z direction is reversed compared to the previous event at tCski = 600 (corresponding to turmskt ~ 300). 




FIG. 2: (Color online) Ux (a) and Uy (b) as functions of time and z 
for Run A. 

FIG. 3: Four times Ux (solid line) and Uy (dashed) from Run A at 
t = 500(csfci)-\ 

In Run A with 512'^ meshpoints there is one particularly 
pronounced event during the time interval 200 < tumis^f < 

400, where Uy reaches an extremum at iuims^f ~ 280, fol- with respect to z are denoted by a prime. 

lowed by an extremum of u'^ a bit later at tu„nskf ~ 300; The occasional extrema in the components of U and its 
see Fig. |4] for their root mean square values. Here, derivatives derivatives are accompanied by strong enhancements in the 
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FIG. 4: Root mean square values of (solid line) and Uy (dashed 
line) for Run A. Note the maxima at tUimski ~ 280 and turmskf « 
260, respectively. 



rms value of the total velocity, Urms, which includes the mean 
flow as well. This fact has been of some significance in pre- 
vious studies of hydromagnetic dynamo action from turbu- 
lent shear flows ifTufTsi 121(1 . because, depending on the value 
of the sound speed, this can lead to numerical difficulties if 
the Mach number exceeds unity during these strong enhance- 
ments of Uinis- These difficulties are here avoided by choos- 
ing a smaller shear parameter Sh, regulated by the input pa- 
rameter S. 

The effect of increasing S is demonstrated in Fig. |5] which 
shows the rms values of the large-scale velocities for four runs 
where S is varied while the other parameters are kept constant. 
The amount of shear is here quantified by the value of Sh, 
which is based on the Uj-ms value from a run without shear and 
thus effectively quantifies the strength of the random forcing. 
These runs are denoted by the letters B to E, with the strength 
of the shear increasing from Sh = —0.08 in Run B to Sh = 
—0.33 in Run E. The bottom panel of Fig.|5]shows that J/ims 
increases almost in proportion to the shear for — Sh > 0.25. 
The flow in the large Sh runs is also highly fluctuating during 
periods of vigorous vorticity generation, see Fig.|6]for a space- 
time diagram of the large-scale velocities from Run E. Even in 
the lowest shear run (Run B with Sh « —0.08), which is very 
similar to the non-shearing case during most of its evolution, a 
weak large-scale pattern is discernible at times, see times after 
tunnskf > 600 in Fig.|2l 



IV. INTERPRETATION 

In order to shed some light on the mechanism responsible 
for the generation of large-scale vorticity, we consider mean- 
field equations ifTil I22I1 . Adopting averages over the {x,y) 
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FIG. 5: Root mean square values of Ux (a) and Uy (b), and ?7rms 
(c) for Runs B to E with different shear as indicated in the legend 
in panel (a). The Reynolds number based on the Urms from a non- 
shearing run is ~ 24. 



plane, denoted here by an overbar, we have 



dU 

dt 



-SU^y + T + vU , 



(6) 



where — —u ■ Vm is a term that results from the nonli- 
nearity of the Navier-Stokes equations, and primes denote a 
z derivative. Note that we have assumed solenoidality, i.e. 
V ■ J7 = Us^s — 0, so C/3 = const = by a suitable choice of 
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FIG. 6: (Color online) Ux (a) and Uy (b) as functions of time and z 
for Run E with Sh ^ -0.33. 
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FIG. 7: (Color online) Ux (a) and Uy (b) as functions of time and z 
for Run B with Sh ^ -0.08. 



the initial condition. Thus, only the x and y components of U 
are non- vanishing. Therefore, U ■ "VU = 0. Furthermore, the 
pressure gradient term does not enter in Eq. because any 
horizontally averaged gradient term can only have a z com- 
ponent. Using mean field theory ifTill . ^ can be expressed in 
terms of derivatives of the mean flow. In the present case of 
one-dimensional mean fields this relationship reduces to 



(7) 



where Vij is the eddy viscosity tensor. We also assume in- 
compressibility of the small-scale velocity field, V • it = 0, 
which is a good approximation for small Mach numbers, and 
recall that horizontal averages depend only on z, i.e. the j = 3 
coordinate. Therefore we have 



(8) 



Here we have denoted the two relevant components of the 
Reynolds stress tensor by TZi = uiu^, where i — 1,2 refer 



to the X and y directions and U3 is the z component of the 
velocity fluctuation. Integrating Eq. Q over z, we have 



Tli + i^ijUj — const. 



(9) 



Given that TZ and J7 can be obtained from the simulations, 
we can then find all four components of I'ij by considering 
moment equations of the form 



v., 



0, 



(10) 



where we have introduced the correlation matrix Mjk = 

{u'jU'i^}. We have also assumed that is independent of 
z, and that, owing to periodic boundary conditions, the mean 
flow and its z derivatives have zero volume average, i.e. 
(C/j) = for any i. The components of i^ij can then be written 
as 



Vt2 



= -M" 




(11) 



where z = 1 or 2. 

It turns out that the components of the correlations (TZiUj) 
are small compared with those of Mj^- This makes the eval- 
uation of the components of i'ij using Eq. (fTTI ) ill-behaved 
(see Fig. [8]l. This procedure does, however, yield reasonable 
results for the 1^21 components: 1^21 is highly fluctuating, but 
with an average of the order of roughly half of the reference 
value i>to = ^Urmsk^^, whereas 1^22 is positive and between 
one and two times vto in the quiescent phases of the simu- 
lation and peaking at roughly 5fto when the vorticity peaks. 
The definition of j/to is analogous to a corresponding reference 
value for the magnetic diffusivity (l^, but it is not clear that 
ft should be exactly equal to uto in any limit. Instead, accord- 
ing to the first order smoothing approximation, ut = 0.4fto 
iflo l, and hence the magnetic Prandtl number was expected to 
be 0.4. ^ 

If both components of {TZiU^) for /c = 1 and 2 were ex- 
actly zero, we could calculate 1^12/1^11 in terms of the ratios 



V12 
1^11 



Mifc 

M2fc 



(12) 



for k = 1 and 2. Yet another possibility is to take the geomet- 
ric mean of the two expressions, so 



'^12 



In M12 
I21 M22 



1/2 



1/2 



(13) 



where we have used the fact that M21 = M12. The results 
shown in Fig. |9] indicate that the two ratios in Eq. (fT2] i give 
consistently negative values, although their moduli are differ- 
ent. Assuming that vn is positive, which is reasonable, this 
result suggests that a negative V12 is present in the system with 
a modulus that is between 0.2 and 0.4 times the vn compo- 
nent. 

Finally, a completely different approach for obtaining esti- 
mates between the components of Vij is to use the resulting 
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FIG. 8: Components of the eddy viscosity tensor as obtained from 
Eq. jllb normalized by uto = |itrmsfcj"^ for Run A. Note that the 
average value of 1^21 is negative (see the dashed line in the third panel 

for tUrmski > 100). 



mean-field equations, Eq. (|6]l, and apply them to a hypotheti- 
cal steady state. These equations are linear, which is a conse- 
quence of assuming the components of to be constant. In 
that case we can Fourier transform and obtain the two equa- 
tions 

(i^ + i^n)k^Ui + ui2k'U2 = 0, (14) 

{S + V2lP)Ul + {v + V22)k''U2=Q, (15) 

where C/i and U2 are the Fourier amplitudes of the x and y 
components of the mean flow. Since these equations are lin- 
ear, they cannot describe nonlinear saturation of a mean-field 
vorticity dynamo instability. However, it is plausible that the 
assumption of constancy of the components Vij breaks down 
when the resulting mean vorticity has become large enough. 
The resulting modifications of Vij may then explain satura- 
tion. 

Equations (fT4l i and (fTSI ) show that a necessary condition for 
the mean-vorticity dynamo to be excited is that the product 
V12S is positive. This is indeed the case; in our case both 
V12 and S are negative. A sufficient condition for the mean- 
vorticity dynamo to be excited is that the parameter 



D = [ui2{S/k' + V21) + e'] /i^T > 1 
where — v + vt with 

vt = + 1/22), 



V22)- 



(16) 



(17) 



The parameter D plays the role of a mean-vorticity dynamo 
number. The assumption of a steady state in Eqs. ( fT4l i and (flST l 
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FIG. 9: Scatter plots of M12/M22 (a), M11/M21 (b), M11/M22 (c) 
for Run A. 



impHes that D = \. Note that Eqs. (fT4l) and ( fTSI ) yield 



V + V22 



V + vyx U2 S/k^ + V21 



(18) 



This allows us to calculate v + V22 in terms of Urms/fcf, pro- 
vided ^21 is negligible or known: 



V + V22 _ Ui 



Sh 



ki 



V2lkl 



(19) 



The amplitudes Ui and U2 for Run A are shown in Fig. [TOl 
Putting in numbers, U1/U2 = 0.23, kf/ki — 5, we obtain 



= 1.15 - 0.23^^, 

^rms/ "'f ^rms 



(20) 
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FIG. 10: (a): Ua: (solid line), 4 C/^ (dotted) and Uy (dashed) as func- 
tions of time for Run A. (b): scatter plot of Ux versus Uy for the 
same run. The dashed line shows a linear fit to the data. 



so the uncertainty in 1^21 enters only weakly. Note, however, 
that //22 is more than three times larger than the original esti- 
mate of h'tQ. 



V. EDDY VISCOSITY FROM THE IMPOSED SHEAR 

We have so far only looked at the components of the 
Reynolds stress tensor that enter the horizontally averaged 
equations. However, there is at least one other component 
that does not enter Eq. but that can also be used to deter- 
mine the eddy viscosity (see, e.g. ||24ll ). This component is not 

driven by the derivatives of U, but by the imposed shear flow, 
g 

VxUy, itself. Indeed, one expects that this imposed shear 
leads to an xy stress 

/ 5 5\ 

uiPu^ ^ -i^t [yxUy +VyUxj ^ -vtS. (21) 

This is indeed the case; see Fig.[TT] It turns out that the i^t de- 
termined in this way is rather similar to the value of 1^22 esti- 
mated from Eq. ( l20l l. Again, there is no good reason that these 
values are the same, because the eddy viscosity obtained from 
Eq. (l2Tl i belongs to a different component of the full rank-4 
eddy viscosity tensor and is not part of the rank-2 tensor con- 
sidered above. 
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FIG. 11: Turbulent viscosity for Run A, as obtained from the 
Reynolds stress component Rxy, divided by the estimate Pto- 



VI. CONCLUSIONS 

The present work has demonstrated quite clearly that in 
non-helical shear-flow turbulence a large-scale flow pattern 
emerges spontaneously. In the present case, where in the ab- 
sence of shear the turbulence saturates at a Mach number of 
order 0.01, the large-scale flow becomes exceedingly strong 
and saturates at a Mach number of 0.1-0.2. This behavior is 
seen both at smaU and at the largest Reynolds numbers consid- 
ered here (Re = 100, based on the inverse forcing wavenum- 
ber). 

The flow pattern can be particularly well pronounced at cer- 
tain times and shows a long wavelength variation in the di- 
rection perpendicular to the plane of the shear flow (here the 
z direction). For negative shear, the x and y components of 
the shear flow are in phase in a way that is compatible with 
an interpretation in terms of a large-scale vorticity dynamo, 
as explored first by Elperin, Kleeorin, and Rogachevskii jlill . 
This means that the large-scale flow is driven by an anisotropic 
eddy viscosity tensor. Particularly important is its xy compo- 
nent, i^xy, which describes the production of a cross-stream 
large-scale flow component, Ux{z,t), from a z-variation of 
the streamwise large-scale flow, Uy{z, t). The mean-vorticity 
"dynamo cycle" is completed by a suitable action of the shear 
itself, which produces a streamwise large-scale flow compo- 
nent, Uy, from the cross-stream component, Ux, by the term 

-U-VU'^. 

The mean-vorticity dynamo cycle can only work if the sign 

5 

of Vxy is the same as that of the shear, VxUy The present 
investigations suggest that this is indeed the case. However, it 
is desirable to verify the sign of Vxy using a test-flow method 
analogously to the test-field method used in magnetohydro- 
dynamics. Some care in using the correlation method is in 
order, because there are examples in magnetohydrodynamics 
where the correlation method give incorrect values for some 
components of the magnetic diffusion tensor, although other 
components were correct |25]. For example, when we apply 
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a method analogous to that in Eq. (fTsT i to the magnetic field 
of a simulation of shear flow turbulence (see, e.g., Figs. 7 or 
8 of Ref. lO); the components of the magnetic field scat- 
ter almost isotropically about the origin. This is compatible 
with an interpretation in terms of an incoherent alpha-shear 
effect lfl2l. I26I1 . On the other hand, there is still a weak cor- 
relation with a negative slope. This would suggest that the 
shear-current dynamo might also be at work, even though the 
test-field method indicates that this should not be the case. 

Clearly, the reality of the large-scale flow found in simula- 
tions is more complicated than what is suggested by the sim- 
ple mean-vorticity dynamo problem. Firstly, in contrast to the 
magnetic dynamo no kinematic stage can be distinguished, i.e. 
the large-scale patterns are visible only after they are already 
of dynamical importance. Secondly, the mean flow can re- 
verse sign in random intervals which is not anticipated from 
the linear mean-vorticity dynamo model with anisotropic eddy 
viscosity, where self-excited solutions would always be non- 
oscillatory. Another question that needs to be addressed in 
future work is the saturation level of the large-scale flow, its 
relation to the saturation level of the small-scale flow, and a 
possible dependence on the Mach number. 

For more realistic applications it will be important to get 
information about the full eddy viscosity, which is a rank-4 
tensor [8J. In the present work, where the averages are only 



one-dimensional, the eddy viscosity reduces to a rank-2 ten- 
sor Finally, for astrophysical applications it should be pointed 
out that the gas in many shear flows is ionized and electri- 
cally conducting, giving rise to efficient dynamo action. The 
resulting mean Lorentz force from the small-scale magnetic 
field modifies the eddy viscosity in a way that suppresses the 
mean-vorticity dynamo. Details of this need to be investigated 
further Another effect that can suppress the mean-vorticity 
dynamo is rotation ifisll . This can be understood from the dis- 
persion relation in that the addition of rotation leads, among 
other terms, to a —411^ term inside the squared brackets of 
Eq. ( fTSl l that always suppresses the mean-vorticity dynamo. 
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